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ABSTRACT 

Modeling the radiation generated by accreting matter is an important step towards realistic simu- 
lations of black hole accretion disks, especially at high accretion rates. To this end, we have recently 
added radiation transport to the existing general relativistic magnetohydrodynamic code, Cosmos++. 
However, before attempting to model radiative accretion disks, we have tested the new code using 
a series of shock tube and Bondi (spherical inflow) problems. The four radiative shock tube tests, 
first presented by Farris et al. (2008), have known analytic solutions, allowing us to calculate errors 
and convergence rates for our code. The Bondi problem only has an analytic solution when radia- 
tive processes are ignored, but is pertinent because it is closer to the physics we ultimately want to 
study. In our simulations, we include Thomson scattering and thermal bremsstrahlung in the opacity, 
focusing exclusively on the super-Eddington regime. Unlike accretion onto bodies with solid surfaces, 
' super-Eddington accretion onto black holes does not produce super-Eddington luminosity. In our 

examples, despite accreting at up to 300 times the Eddington rate, our measured luminosity is always 
q ' several orders of magnitude below Eddington. 

Subject headings: accretion, accretion disks — black hole physics — magnetohydrodynamics (MHD) 
c/j ■ — methods: numerical — radiative transfer 

^ ■ 1. INTRODUCTION 

After more than a decade of successes in the area of general relativistic ma gnetohydrodynamic (GRMHD) numerical 
OO ' simulations of black hole accretion flows (see lAbramowicz fc Fragile! 1201 ll for a review), there is now considerable 
, understanding of accretion in this context. However, one context in which our understanding is still quite limited, is the 
■ inner parts of luminous accretion flows. Here radiation pressure is expected to dominate over gas pressure in supporting 
' the flow against the gravitational field of the black hole. Because of photon diffusion, such radiation-pressure-dominated 
\ regions can be even more fragile to the development of inhomogeneous density structures than gas-pressure-dominated 
f^S • ones, and because of the proximity to the black hole, relativistic effects are of utmost importance. For these reasons, 
£N| ■ numerical study of this and similar problems demands the development of general relati vistic radiation MHD codes. 

Although radiation hydrodynamics codes have been around for more than 4 decades (e.g. lChristvll 966: Colgate fc White! 
kJ ! [1961 ICox etalWwM ].' very few multidimensional, fully (general ) relativistic treatments exist ev e n tod a y. A few note- 
worth y examples closely related to our work were presented bv iDe Villiersl (|2008l ): iFarris et al.l (|2008[ ): iZanotti et al.l 
(|20T1 . Parallel progress has also been made in developing radiation hydrodynami cs codes to solve the problem of 
core-collapse supernovae ([Muller et al.l [20Toi : iShibata et al.l I20TH iLentz et alll2012T) . There the radiation is in the 
form of neutrinos, but the processes (and code challenges) are very similar. The paucity of codes can be understood 
when one realizes that most radiation schemes become quite difficult to implement in multidimensions, given the large 
number of phase space degrees of freedom that need to be tracked for the radiation field. The sheer computational 
demand of multidimensional radiation MHD has, until recently, been beyond the reach of even the largest scientific 
computers. 

The potential for discovery from such simulations, though, is significant. The addition of radiative phenomena into 
general relativistic fluid dynamic simulations greatly expands the possibl e parameter space that ca n be explored, with 
interesting phenomena such as radiation- press ure-dominated s lim disks (jAbramowicz et aDll988f ). secular instability 
([Lightman fc Eardlevl lT974). photon bubbles (|Begelmanl [200 D . and radiation driven outflows, all topics that could 
potentially be investigated. Therefore, in this paper, we take the first, simplest step, c onsidering spherical accretion 
onto a non-rotating black hole. This case has been considered many times previously ([Shapiro! 119731 : iSchmid-B urgk 
ll97aiMeienll979l:lSoffell982l:lVitellolll984l:lZampieri et aLlll996l ). but never, to our knowledge, with an explicit, multi- 
dimensional, Eulerian, general relativistic radiation MHD code such as the one implemented here. For this reason, we 
feel it is worthwhile for us to present these results. 

Although our implementation of GR radia tion MHD is not unique - in fact, it is very similar to at least two previous 
implementations described in the literature ([Farris et a l. 2008; Zan otti et al.ll20lTI ) - we, nevertheless, begin our paper 
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in fJ2]with a detailed discussion of our method. In $3] we report on a series of one-dimensional shock tube tests meant 
to validate our code. In $4] we arrive at the main results of this paper - one-dimensional spherical accretion onto a 
black hole including radiation. We conclude in <J5] Most of the equations in this work are written in units where 
GM = c = 1, although in a few places we leave in factors of c for clarity. 

2. NUMERICAL METHOD 
2.1. GR MHD 

We begin with a review of general relativistic MHD before discussing the addition of radiation. This will make it 
easier for us to point out the advantages of our implementation. We start with the three basic laws of GRMHD: the 
conservation of the stress-energy tensor 



a/3 



the continuity equation 

and the homogenous Maxwell's equation 



T 







o , 



In the MHD limit, the stress-energy tensor can be expressed as 

T° p = (ph + 2P niag )u a u fj + (P ga 



P mag )s^ - b a b? 



(1) 
(2) 
(3) 

(4) 



where p is the rest mass density, h = 1 + e + P gas /p is the specific enthalpy, e is the specific internal energy density, 
P gas is the gas pressure, P mag is the magnetic pressure, u a = g a/3 up is the fluid 4-velocity, g a p is the 4-metric, g is the 
4-metric determinant, F a/3 is the Faraday tensor, and *F a ^ — e a P lS is its dual. For this work we use a T-law 

equation of state (EOS), 

P=(T-l)pe (5) 

where T (without subscripts or superscripts) is the adiabatic index. 

Expanding the MHD equations and reorganizing in terms of new variables, we arrive at the following set of coupled 
equations 



d t D 
d t £ + d, 
d t Sj + d, 
d t B j + d l (B j V 



d l {DV l ) = Q , 



B*V J ) = 



(6) 
(7) 
(8) 
(9) 



where D = Wp is the generalized fluid density, W = yj^gu = ^/^gj/a, 7 = l/Vl — v 2 , a = l/y— 3 00 , v 2 = ViV 1 , 
v % is the flow velocity relative to the normal observer, V % — u l /u° is the transport velocity, £ = —^J—gT® is the 
total energy density, and Sj — yJ—gT® is the covariant momentum density. With indices, T indicates the geometric 
connection coefficients of the metric. Along with the evolution equations ((SJ-©, the time component of the Maxwell 
equation enforces the divergence-free constraint djB J = 0. 

There are multiple representations of the magnetic field in our equations, which we should explain: b a is the magnetic 
field measured by an observer comoving with the fluid, which can be defined in terms of the dual of the Faraday tensor 
b a = up*F a P , and B^ = ^J—gBi is the boosted magnetic field 3-vector. The magnetic field B % — *F m is related to 
the comoving field by 



B' 



.Ohi 



mO 



(10) 



We also need the magnetic pressure, which is defined as P mag = b 2 /2 = b a b a /2. Note that unlike previous versions of 



Cosmos++, here we have absorbed the factor of \/4ir into the definition of the magnetic fields. 

2.2. HRSC 

The current work uses a high- resolution shock capturing (HRSC) scheme to solve the GRMHD equations ©-(HJ). 
This method is ne w to Cosmos++, although it shares many elements with our earlier non-oscillatory central difference 
(NOCD) method JA nninos fc Fragile! [2003D . We begin by noting that the MHD equations are all in the form of a 
conservation equation 

a t u(P) + a,F i (p) = s(P) (11) 



where 



U(P) = 



/ D 

£ 



F ,; (P) 




S(P) = 









\ 



(12) 
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are the arrays representing the conserved quantities, fluxes, and source terms, respectively. By integrating both sides 
of equation (JTTJ) with respect to volume and applying Gauss' Theorem, we can rewrite the conservation equation in 
the form . r r 

/ dtUdV = - i F l dA l + / SdV , (13) 
Jv Js Jv 

where we have dropped the explicit dependence on the primitive variables. We can then discretize this equation using 
a finite volume representation as . 

u n+1 = u n - — ( FM + Ats • ( 14 ) 

faces 

This approach requires at least a 2nd order time integration scheme for stability. The Cosmos++ code has a number 
of time integration options, including: 2nd order Euler, 2nd order Runge-Kutta, 2nd order Crank-Nicholson, and 3rd 
order Euler. The present work uses the 2nd order Runge-Kutta scheme (which is also sometimes referred to as the 
"midpoint" or "leapfrog" method). First a half-time step, At/2, is taken to project the conserved variables U" forward 
to n + 1/2. From these, a new set of primitives p n + 1 / 2 can be computed. These intermediate primitives are then used 
in calculating the F* and S needed in equation (|T3|). 

The F* are determined using an approximate Riemann solver. We have options for either the HLL or Lax-Friedrich 
method. The HLL scheme reconstructs the fluxes as 

F = . (15) 

Cmax T C m i n 

A slope- limited linear or parabolic extrapolation gives Pr and Pl, the primitive variables at the right- and left-hand 
side of each zone interface. From Pr and Pi,, we calculate the right- and left-hand conserved quantities (Ur and Ul), 
the fluxes Fr = F(P^) and F^ = F(P^), and the maximum right- and left-going waves speeds, c± t R and c± t L- The 
bounding wave speeds are then c max = max(0, c + .r, c + x) and c m ; n = — min(0, c^u, c_.l). By setting c max — c m i n , 
the HLL flux can be reduced to the local Lax-Friedrichs flux. 

For the parabolic interpolations of P# and P^, we use the piecewise-parabolic method (PPM) of lColella fc Woodward! 
()1984l) . As noted by those authors, this method can, on occasion, produce unwanted oscillations, especially behind 
stationary s hocks, such as the ones to be treated in $3l To combat these oscillations, we added the flattening procedure 
described in lColella fc Woodwardl (|1984D . 

As noted above, the conserved variables U, the fluxes F*, and the source terms S are all functions of the set of 
primitive variables 

/ P \ 

(16) 




One difficulty in relativistic MHD is the so-called "inversion problem," that is going from updated conserved variables 
to updated primitive variables. Unlike Newtonian MHD, there is not a set of analytically-solvable algebraic expressions 
for this inversion. Instead, one must use a numerical procedure, such as the Newton- Raphson method, to solve one or 
more of the inversion equations. In our code we have implemented the 2D, IDw, and ID V 2 methods of iNoble et al.l 
(2006). Our code defaults to the 2D method, but will fall back to the other methods in succession if the 2D method 
fails to converge. 

2.3. Constrained Transport 

The finite- volume discret i zation of the induction equation presented in £12.21 can be treated using the constrained 
transport schemes of iTothl (|2000| ). However, these methods have certain inadequacies we wish to avoid. Most im- 
portantly for us, they are not easily extendible to adaptive mesh refinement. Another significant shortcoming is 
the rather large stencil they require. The methods can also lead to unphysical behavior for certain types of flows 
(IGardiner fc Stond[2008h . 

Instead, we depart from the volume-averaged representation of all of the fields to use a staggered representation of 
the magnetic fields, with the primary representation of the fields being face-centered. We then integrate the induction 
equation (j9]) over surfaces rather than the volume of the cell. For example, integrating equation ([9]) over the x\-, xi~, 
and X3-faces located at (i — 1/2, j, k), — 1/2, k), and (i,j,k — 1/2) respectively gives, after use of Stokes Law, 
(|Stone fc Gardinedl2009l) 

1° k-l/2.j,k-y& H-X/2,j,k ~ -faT (< t, 3)i-l/2.i + l/2.k ~ l®3j i _ 1/2 .-i-l/2.fe I 1 ' > 



St 

Sx 3 



I e \«+l/2 / m \ri+l/2 

y^)i-l/2,j,k+l/2 V & 2)i-l/2,j,k-l/2 



5t 



(n2\n+l —(K 2 \ n "' (e \ n + 1 / 2 le \"+l/ 2 

y' 6 >i,j-i/2,k - ytt >i,j-i/2,k - [v®iJij-i/ 2 ,fc+i/2 ~ \®Vi,j-i/2,k-i/2\ y 1 *) 



St 
Sx\ 



i e \ ft+l/2 / n, xn+1/2 

V 6 3j l+1/2 ,,-i/2.fe ~ y^lt-l 



i+l/2j-l/2,fc y" A h-l/2,j~l/2,k 
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St 



,K : M"+' _ , R-'t ,» / » \n+l/2 

V° ^,j,fc-l/2 — ^ >i,j,k-l/2 ~ fa^ |_l 6 2j l+1/2 j ife _i /2 



/ «. xn+1/2 
l <S 2j l _ 1/2j:fc _ 1/2 



(19) 



St 
5x 2 



\ 01 )i,j+l/2,k-l/2 



I e \«+l/ 2 
\ &1 )i,j~l/2,k-l/2 



where, for example, (A)ij-i/2,fc-i/2 is the a; i -component of the electric field (or emf) Si — — eyfeV^ B k centered on 
the appropriate cell edge. The best method we have found so far for constructing the edge-centered electric fields is 
to simply average the surrounding face-centered fluxes recovered from the Riemann solver. For example, 

1 



0*0 



i-l/2,j-l/2,fc 



(<%) 



i-l/2J,k 



/2,j-l,k 



0*3 



i,j-l/2,k 



+ (<%) 



i-1 ,j-l/2,k 



(20) 



This procedure preserves the following mathematical representation of the divergence to round-off error in each zone: 



1 



-{B 2 A 2 )ij-i/2,k + (S 3 A 3 ) i j :fe+1 / 2 - {B Az)i t j t k-l/2~\ 



(21) 



In this scheme we maintain staggered forms of both the primitive and conserved magnetic fields, which are related 
by the metric determinant at their respective faces, e.g. B 1 = {^/—g)i~i/2j,kB 1 . The staggered primitive field is used 
to overwrite the appropriate component of the extrapolated primitive field used in the flux reconstruction at each face. 
They are also used when we require cell-centered values of the primitive magnetic fields, such as for calculating the 
magnetic pressure. In this case we use the volume-averaged fields 

[{B )i-l/2,j,k + (B 1 ) i+ l/2,j,k\ 



(B )i,j,k ■ 
■ I'' 2 ),.,.*■ : 

(B )ij,k ■ 



[{B )jj-1/2,& + {B 2 ) i j + 1 / 2 ^} 



(22) 
(23) 
(24) 



[(B 3 )ij. k -l/2 + {B )i,j,k+l/2] ■ 

2.4. Radiation 

We now describe the addition of the radiation field. Similar to the fluid and magnetic field, the energy and momentum 
of the radiation field are represented by a stress-energy tensor 

R af3 = J I u n a n l3 di'dn , (25) 

where v is the photon frequency, I v = I(x a ;n l ,v) is the specific intensity of radiation at position x a moving in the 
direction n a = p a /hpv, p a is the photon 4-momentum, hp is Planck's constant, and dfl is the differential solid angle 
around the direction of propagation. The quantities v, I u , and d£l are measured in the local Lorentz frame of an 
observer with 4- velocity u" &d ^ . 

GR radiation MHD, then, starts from the same conservation of stress-energy tensor as GRMHD, but now split into 
MHD and radiation components: 

(T afs + R af) ) = . 



It is more convenient for our purposes to rewrite equation (|26[) in two parts: 



and 



(T Q/3 ) = G a 



[R^). = -G a , 



(26) 

(27) 
(28) 



where G a is the radiation 4-force density coupling the fluid and radiation field. 

2.4.1. Radiation Moments 

Similar to fluid dynamics, we can proceed by considering various moments of the radiation equations. Starting with 
the zeroth radiation moment, p 

E= I v dvdtt (29) 

is the comoving radiation energy density, . 



l v ridvd£l 

is the comoving radiation flux (three first radiation moments), and 



I u n l n 3 dvdfl 



(30) 



(31) 



is the (symmetric) comoving radiation pressure tensor (six second radiation moments). In the comoving frame, the 
radiation stress tensor is then 



R &0 = 



( E 


F 1 


F 2 


F 3 


F i 


<pii 


<pl2 


<pl3 


F 2 


p21 


j,22 


j,23 


\ F i 


p31 


p32 


p33 



(32) 
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In our current implementation, we assume the radiation pressure, -P ra d, is isotropic, such that V 1 ^ = P ra d = S^E/3. 
This is equivalent to keeping the first two radiative moment equations, and adopting an Eddington factor of 1/3 to 
close the set. Although the radiation pressure is assumed to be isotropic, we emphasize that by allowing a small, 
nonzero radiation flux (F a <C E), we do allow for some degree of anisotropy in the radiation field. In covariant form, 
R a P can now be written as 



R a P = Eu a u p + F a U P + F 



Pmdh 



a0 



(33) 



where we have introduced the tensor h a @ = g al3 + m q u' 3 , which projects any other tensor into the space orthogonal to 
u a , such that h a ^u a — 0, and the radiation flux 4- vector 

F a = h a p J I u n' 3 dvdVL . (34) 

With this definition, the flux satisfies F a u a = 0. 



2.4.2. Radiation Four-Force Density 

We return now to consider the form of the radiation 4- force density from equations (|27p and (|2"5|) . In the comoving 
frame 

G & = J {xuh - ri„) n & dvoin , (35) 

where Xv = xl+xl is the total opacity ( "a" and "s" stand for absorption and scattering, respectively) and r\ v = rft+rfl is 
the total emissivity. In terms of the usual cooling function A = J r\ v dvd£l and assuming a mean (frequency-independent) 
opacity of the form \ = pn, we can write in the comoving frame 



G° 



{Xvlv — f]v) dvdQ = pn a E — A , 



and 



G l = / {Xvlv ~ lb,) n l disdn = pnF 1 = p(n a + k^F* 



(36) 



(37) 



Equation (1371) assumes that photons are emitted isotropically in the comoving frame, so that the net momentum they 
remove from the gas is zero. 

Noting that u a = (1,0,0,0), then, from the normalization of the radiation flux 4-vector F a u a — 0, we have 
F a — (0,F l ). Thus, in the comoving frame 



G & = (pn a E - A) u & + P kF & . 



(38) 



Since this expression is covariant, it must hold in any frame; we, therefore, drop the " in all subsequent references to 
G a . 

For equilibrium blackbody radiation, A = pK a a^Tg as , where T gas = mP gas /kBP is the ideal gas temperature of the 
fluid, aji — 8ir 5 kg / (WhpC 3 ) is the radiation constant, ks is Boltzmann's constant, and m = /xmn is the mean mass of 
ions in the gas. If the radiation and gas are in local thermodynamic equilibrium, then E — a^Tg as and the first two 
terms of equation ()38[) drop out, though we emphasize that our scheme does not require this. 



2.4.3. GR Radiation MHD 

Again expanding and reorganizing the conservation of stress-energy in terms of new variables, we get four new 
conservation equations for the radiation 



d t TZ + ft 
d t TZ, + a 



9 Rp r 0a 



■tn,j t ui vv -.9 Rj) = V^g Rp r^ Q 



-g G 

-a G i > 



(39) 
(40) 



where 1Z = ^—gR® is the conserved radiation energy density and TZj = ^—gRj is the conserved radiation momentum 
density. The addition of the radiation fields also modifies the MHD energy-momentum conservation equations. The 
new equations are 



d t £ + ft (- V=5 n) = -V=g t$ r<L 

dtSj + ft {J—g Tj) = y/=g Tp rf Q + 



The full set of conserved variables, fluxes, and source terms is now 

F*(P) = 



I £ \ 



U(P) = 



K 



DV l 
y/=9T\ 

V=g~Rh 
^gR) 



S(P) = 



G , 
-9 Gj , 





~g t% rft, 
-g t% r 



(41) 
(42) 
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jot 



P 



V~-g R p r oa 

V sTg-R^T 



P ja 



_ \ 

-g G 
~gG 3 

~g Go 
~gG, ) 



(43) 
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The expanded set of primitive variables is now 



I P \ 

V J 
D 1 
E 

V FJ J 



(44) 



Fortunately the inversion of the radiation variables can be written as a simple set of algebraic expressions (utilizing 
-frad = and can be solved independently of the MHD primitive fields: 



E = 3 



F" 



TZ° + 2u°u a TZ a \ 
s/^gg 00 - 2Wu° J 

(WE + u a K a ) 



W 



-Eu J 



9 



oj E 



3u° 



(45) 
(46) 



(47) 



Inside the code we check that the radiation flux satisfies the physical limit J ' F^F^ < E. 

The advantage of this approach is that, by writing the radiation equations in conservative form, we can take full 
advantage of the HRSC machinery described in §2.21 to solve the full set of GR radiation MHD equations. Thus, there 
is relatively little code development beyond adding new variables, flux terms, and source terms. 



2.5. Wave Speed Calculation 

The HLL and Lax-Friedrich approximate Riemann solvers only require the maximum and minimum wave speeds, 
as opposed to the full eigenvectors of the characteristic matrix as would be necessary for a Roe-type scheme. These 
wave speeds are also required to fix the time step via the Courant condition. The relevant speed is the phase speed 
uj/k of the wave, and we treat each coordinate direction independently. For signals propagating in the x\ coordinate 
direction, the corresponding eigenvector is k a = (— u>, k\, 0, 0), and the wave speed is iojk\. 

To find the necess ary wave speeds in the grid frame, we s tart with the following approximate dispersion relation in 
the comoving frame ([Gammie et al.ll2003t iFarris et al.ll2008f ) 



uj: 



v T k, 



i 

cm 



where 



v\ = b 2 /(ph 
relations, uj rrn 



i4 



1/3 



v 



1 



(48) 



(49) 



b 2 ) is the Alfven speed, and c 2 = TP gas /ph is the sound speed. We then substitute the following 
= — k n u a , k 2 ,„ = K n K a , and K a = (g a p + u a up)kP , into equation (j4"8)) to get the following quadratic 



equation for the desired wave speed w/ki 



,00 



^ +2 



t i V 1 

vW 1 



Wave speeds in the and X3 directions are found analogously. 



„oi 



V 1 



vW 1 + 



(50) 



3. RADIATION SHOCK TUBE TESTS 

Unfortunately, there are very few good test problems for r elativistic radiation MHD codes at this time. Among 
the few are the four radiative shock tube tests introduced in IFarris et al.1 (|2008| ) , which we now use to validate our 
code. Each test includes a nonlinear radiation-hydrodynamic wave: case 1 is a nonrelativistic strong shock; case 2 is 
a mildly relativistic strong shock; case 3 is a highly relativistic wave; and case 4 is a radiation-pressure-dominated, 
mildly relati vistic wave. 

Similar to lZanotti et al.l ()2011[ ). and unlike IFarris et all (|2008j), we initiate these problems as traditional shock tubes 
with two states ("Left" and "Right"), initially separated by an imaginary partition, instead of starting from the 
analytic solution. At t = the partition is removed and the gas and radiation are allowed to evolve until a steady state 
is reached. The initial parameters for the four tests are presented in Table[T] The initial values for the flu xes, not given 
in Tab led] are se t to F x = \Q~ 2 E. The scattering opacity k s is set to zero in all these tests. Also like lZanotti et all 
(|201lD and unlike IFarris et a l. (2008), we do not boost the fluid velocities, instead presenting our results in the frame 
of the shock. This is important, as quasi-stationary shoc ks present a particular diffic ulty for the PPM scheme we are 
using. In fact, without using the flattening procedure of iColella fc W oodward (1984), we find that the two tests that 
exhibit discontinuities (cases 1 and 2), suffer from post-shock oscillations with amplitudes of a few percent. With the 
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Table 1 

Shock Tube Parameters 



Case 


r 




PL 


Pl 




u x L 


E L 


PR 


Pr 


U R 


Er 


^stop 


1 


5/3 


0.4 


1.0 


3.0 x 10- 


-5 


0.0015 


1.0 x 10" 8 


2.4 


1.61 x 10~ 4 


6.25 x 10~ 3 


2.51 x 10" 7 


4000 


2 


5/3 


0.2 


1.0 


4.0 x 10" 


3 


0.25 


2.0 x 10" 5 


3.11 


0.04512 


0.0804 


3.46 x 10" 3 


3000 


3 


2 


0.3 


1.0 


60.0 




10.0 


2.0 


8.0 


2.34 x 10 3 


1.25 


1.14 x 10 3 


100 


4 


5/3 


0.08 


1.0 


6.0 x 10" 


3 


0.69 


0.18 


3.65 


3.59 x 10~ 2 


0.189 


1.30 


500 



Table 2 

L-l Norm Errors for Case 4 



Grid 


\E(p)\i 


\E(P)\i 


\E(V*)\i 


\E(E)\i 


\E(F*)\i 


1600 
3200 
6400 


1.69 x lO" 6 
4.22 x 10~ 7 
1.06 x 10~ 7 


2.85 x 10" 8 
7.06 x 10" 9 
1.78 x 10" 9 


1.17 x 10" 6 
2.89 x 10" 7 
7.16 x 10" 8 


8.98 x 10" 7 
2.23 x 10" 7 
5.54 x 10" 8 


1.17 x 10" 7 
2.91 x 10" 8 
7.27 x 10" 9 


Convergence 11 


1.99 


1.99 


2.01 


2.01 


2.00 



a Convergence rate between 3200 and 6400 zone data. 



flattening procedure, the oscillations are effectively removed, without creating unnecessary dissipation in the smooth 
test problems. 

Figures [J - [H sh ow the four cases. Although it is possible to calculate semi-analytic solutions for each of these 
(|Farris et al.ll2008t ), we have chosen instead to simply plot results using 800 zone resolution against results using a 
much higher (3200 zone) resolution. In all cases, the results agree quite well at the different resolutions, and with 
previously published results. 

We have used the smooth wave in case 4 to test the convergence rate of our numerical scheme. Table [5] reports 
the L-l norm error (i.e. |S(a)|i = ^2 { Ax\ai — A{\, where a% and Ai are the numerical and semi-analytic solutions, 
respectively) for 1600, 3200, and 6400 zones resolution. The convergence rate of the errors for all variables at all 
resolutions is almost exactly 2, as expected for a smooth flow using our overall scheme. 

4. BONDI INFLOW WITH RADIATION 

As we mentioned in <JIJ the case of optically-thick, spherical accretion onto a non-rotating black hole has been 
considered many times in the past, although usually with strictly one-dimensional codes. Another difference is that 
most of the codes in the past have used implicit integration schemes, whereas we are, for now, attempting to proceed 
using the explicit scheme already in Cosmos++. Finally, instead of assuming local thermodynamic equilibrium between 
the gas and the radiation, we consider two physical cooling processes in the gas: Thomson scattering and thermal 
bremsstrahlung. The first contributes an opacity 

k s = 0.4 cm 2 g- 1 . (51) 

while the second has the form ()Rvbicki fc Lightmanl fl986') 

n a = 1.7 x 10- 25 T K 7/2 p cgs m p 2 cm 2 g" 1 , (52) 

where Xk is the ideal gas temperature of the fluid in Kelvin, p cgs is the density in g/cm 3 , and m p is the mass of a 
proton in g. We assume the gas is fully ionized hydrogen, so the mean mo lecul ar weight ji = 0.5. 

We have chosen to set up cases similar to ones presented in lVitellol ()1984[ ) and lNobili et al.l (|1991[ ). th ough significan t 
differences in our assumptions mean we do not expect to exactly reproduce their results. In the case of lVitellol (|1984h . 
they assume LTE throughout the flow, such that they do not need to solve the radiation transport independent 
of the fluid transport. As we emphasized in £12.4. 2[ we do not enforce LTE in our method. There are even more 
differences between our approach and that of iNobili et al.l (|1991l ). First, they consider a more physical, though also 
more complicated, equation of state. They also account for Compton scattering of the radiation, which we do not 
consider in this work. Finally, they use a constant value of T = 10 4 K for the temperature of the gas at the outer 
boundary of their simulations, whereas we explore values in the range 10 5 -10 7 K. 

We use a logarithmic radial coordinate of the form x\ = 1 + m(r/rs) to cover the expansive spatial range required, 
where rg = 2 is the Schwarzschild radius. All simulations use a resolution of 512 zones. Note that because this problem 
is spherically symmetric, we have performed our tests using a ID, spherical-polar (radial) grid. Nevertheless, we have 
confirmed that our code produces identical results when the problem is run in 2 and 3 dimensions (when using the 
same radial grid spacing). 

To initialize the problems, we first fix the density, p Ql and temperature, T , of the gas at r Q . Once T Q and p are 
fixed, and assuming some relation between T gas and T ra d, we can determine the polytropic index of the gas at r Q from 
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Figure 1. Profiles of p, P, V x , E, and F x at t = 4000 for case 1. Symbols denote data from 800-zonc simulation (sampled to only display 
200 points); solid lines denote data from 3200-zone simulation. 

From this expression, we can easily see that Y = 5/3 for gas-pressure-dominated flows, whereas Y = 4/3 for radiation- 
pressure-dominated ones. We assume the initial value of Y found at r a applies throughout the flow for the duration of 
the simulations. By assuming a polytropic equation of state, F cc /, we can determine the initial temperature profile 
of the gas from 

T sas =T (p/p f- 1 . (54) 

We still n eed to specify the initial profiles of u r and p. For simplicity, we assume that these equal their free-fall values 
u r = —y/2M/r and p = —M/4irr 2 u r at all radii, with the mass accretion rate,_M, now one of our free parameters. 
We explore mass accretion rates in the range 10 < m — M/M^dd < 300, where MEdd is the Eddington mass accretion 
rate. The lower limit is set by the requirement that at least some portion of the flow be optically thick. The upper 
limit is set by the requirement that the "photosphere" be contained within the grid. The photosphere is where the 
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Figure 2. Profiles of p, P, V x , E, and F x at t = 3000 for case 2. Symbols denote data from 800- zone simulation (sampled to only display 
200 points); solid lines denote data from 3200-zone simulation. 



optical depth r becomes 1, with r defined as 

x 

(X a + X s ) ds . (55) 

In practice we approximate the optical depth as r ~ (x a + x s ) r - 

We consider three types of simulations, all with M = 3M Q , r Q — 10 4 rs, and rj = 0.95rs. In the first, the gas is 
purely adiabatic, as we ignore radiation; in the second, we include radiative processes as described above; in the last, 
we also include a radial magnetic field of the form B r = —deA^, where 

A<p — sign(cos 9) cos0y/2P o /l3 o r 2 o (56) 

and p o = 100. Since the field is purely radial and weak, it has no noticeable effect on the dynamics. Nevertheless, we 
wanted to include at least one test using all of the physics discussed in this paper. Future work will rely much more 
heavily on the full capability. 
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X 

Figure 3. Profiles of p, P, V x , E, and F x at t = 100 for case 3. Symbols denote data from 800-zone simulation (sampled to only display 
200 points); solid lines denote data from 3200-zone simulation. 



For the tests with radiation, we initially specify the ratio of radiation to gas temperature, T la a/T gas <C 1. This is 
mainly done so that the radiation energy density E = aT^ ad may start with some reasonable value. The radiation flux 
is similarly set to an arbitrary initial value F r -C E. We confirm that our final results are not sensitive to our choices 
for these parameters. The reason we initialize the radiation temperature to a much lower value than the gas is because, 
if the radiation and gas are nearly in thermal equilibrium, i.e. E s» orT^, then the first two terms of equation Q38p are 
approximately equal, meaning their difference can be arbitrarily small. This makes the corresponding source terms of 
equations p9[) - (|42[) incredibly stiff, and stable evolution with an explicit scheme would require an unacceptably small 
timestep. In future work we plan to explore using an implicit scheme to solve this source term, which will alleviate 
this stability problem. However, this limitation does not significantly affect the simulations presented here. This is 
because the radiative processes we are considering are not very efficient, so the steady-state radiation temperature 
(and pressure) are naturally significantly less than the gas temperatures (and pressures) we consider. 

Table [3] summarizes the key simulation parameters for this section. Each simulation is run long enough for E and F r 
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Figure 4. Profiles of p, P, V x , E, and F x at t = 500 for case 4. 
200 points); solid lines denote data from 3200-zone simulation. 



Symbols denote data from 800-zone simulation (sampled to only display 



to achieve steady-state profiles out to the photosphere. Our main diagnostic in this section is the emitted luminosity, 
which can be recovered directly from the evolved radiation fluxes F % . Specifically, 



L = 



F l n,dA 



(57) 



where A is the surface area encompassing the volume V. The surface is taken to correspond to the photosphere (r = 1). 
We measure the luminosity in units of the Eddington luminosity I = L/L^di where L^d = MsddC 2 = ^kGMcot jrn v 
represents the limit at which outward radiation pressure balances gravity. Above this limit, the radiation pressure is 
sufficient to halt, or even reverse, accretion. We define the radiative efficiency of our flows as ry = l/m. 

Figure [5] shows profiles of six different simulations exploring different values of T D and m. We did not include the 
adiabatic (ElOa) nor radiation + magnetic field (E10T6b) cases in this plot since they are practically indistinguishable 
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Table 3 

Radiative Bondi Simulations 



Simulation 


m 


To (K) 




a 

3 


^stop 


cycles 


I 




E10T5 


10 


10 5 


1.2 x 10" 


7 


7 x 10 3 


7.8 x 10 6 


2.01 x 10- 


6 


E10a b 


10 


10" 






7 x 10 3 


2.7 x 10 5 






E10T6 


10 


10" 


1.2 x 10" 


4 


7 x 10 3 


7.8 x 10 6 


3.72 x 10- 


-6 


E10T6b c 


10 


10 e 


1.2 x 10" 


4 


7 x 10 3 


7.8 x 10 6 


3.72 x 10- 


-6 


E10T7 


10 


10 7 


0.12 




7 x 10 3 


7.7 x 10 6 


1.13 x 10" 


-5 


E30T6 


30 


10 e 


3.9 x 10" 


-5 


7 x 10 3 


7.7 x 10 6 


2.66 x 10- 


5 


E100T6 


100 


10 6 


1.2 x 10- 


-5 


1.6 x 10 4 


1.8 x 10 7 


1.97 x 10- 


-4 


E300T6 


300 


10" 


3.9 x 10- 


-6 


4.3 x 10 4 


4.8 x 10 7 


7.91 x 10" 


-4 



a Measured at the inner radial boundary at t = 0. 
k Adiabatic 

c Radial magnetic field case 



from our reference radiation simulation, E10T6. The fact that the adiabatic and radiation simulations are so similar 
is not surprising since the radiation pressure is so much smaller than the gas pressure, P ra d < 0.025P gas throughout 
the flow; so even with radiation, the flow behaves nearly adiabatically. This reflects the very inefficient nature of the 
radiative processes we are considering. The radiation + magnetic field simulations appears nearly identical, which 
again is expected since the magnetic field is purely radial and weak, so it can not play a dynamical role. 

In all of the cases the inflow is supersonic over the entire radial domain. The accretion radius r a — GM/c^ ^ also 
lies beyond the outer radial boundary for each of these cases. We are, however, able to identify the trapping radius 
r t , where the advection of photons inward becomes faster than their diffusion outward, i.e. |£?V' r | > F r , in all of our 
simulations; this point is marked with a symbol on each of the plots of F r . 

Simulations E10T5, E10T6, and E10T7 illustrate one of the main drawbacks of our current method - the restriction 
to optically thick flows. For simulations with m < 10, the flow is only optically thick in the inner few rg. Once the 
optical depth drops below r 0.1, we begin to see oscillations or noise in the profiles of E and F r . These oscillations 
do not appear to damp away with time, and indicate a fundamental limit to our method. 

Even though the mass accretion rates in all these simulations are highly super-Eddington (m > 10), the luminosities 
are not (I < 7.9 x 10~ 4 ) (see Table[3]). Clearly the low radiative efficiency of this flow 77 = l/fn < 2.6 x 10 -6 indicates 
that not all of the binding energy that is liberated by accretion is able to escape in the form of radiation. Instead, 
much of the energy, in the form of kinetic energy, heat, and radiation, is advected into the black hole. Figure [6] presents 
all of our Bondi results in the l-rh plane. These results are broadly consisten t with earlier studies of optically thick 
spherical accretion (e.g. lVitellolll984t iNobili et al.1ll991t iZampieri etaD 1996). given that different assumptions were 
used in each work. 



5. CONCLUSION 

In the study of black hole accretion, it only makes sense to treat the whole problem (gravitation, gas dynamics, 
magnetic fields, and photons) in a relativistic framework. Therefore, the development of a fully relativistic radiation 
MHD numerical code is an essential step. Along with allowing us to answer some of the very longstanding questions 
concerning black hole accretion, having this capability will also allow us to explore such novel effects as disk self- 
illumination due to light bending near the black hole. There are also astrophysical applications of a relativistic 
radiation MHD code beyond black hole accretion, including core-collapse supernovae, collapsars, and gamma ray burst 
sources, meaning this work will potentially open up new avenues of research. 

This paper demonstrates that we have taken the first step toward that goal. We now have a scheme that is able 
to self-consistently treat the radiation on equal footing with the MHD in a fully general relativistic framework. We 
have sh own th at the method performs well when restricted to the appropriate parameter space, i.e. we require 
\F\ = ■J 1 F^F^ ■C£,t>1, and that the gas and radiation be far from local thermodynamic equilibrium. 

One improvement to our method might be to extend it to treat both optically thick and thin. flows by implementing 
a more general closure relation. Currently, we assume the radiation pressure has the form P y = VS^ — E/iS 1 ^ , i.e. 
we adopt an Eddington factor of 1/3 to close the set of radiation moment equations. Although the simplicity of this 
prescription is desirable, this approximation is only appropriate in the optically thick limit. More general Eddington 
tensors, that rec over the correct asymptotic behavior for both optically thick and optically thin gas, are available (e.g. 
lLevermore]|1984D and could be implemented. 

A related issue that must be addressed is that HRSC schemes fail to treat the optically thick limit properly when 
the photon mean-free path is shorter than the numerical grid spacing (|Lowrie fc Morelll200lD. The simplest fix for this 



is to phase out the relevant terms in the diffusion limit when the mean- free path is short ( Jin fc Levermore1ll996l ). 

Finally (and not surprisingly) , we have found that the fully explicit method described in this paper suffers from severe 
timestep restrictions. Even the relatively modest 512 zone, ID spherical simulations in S|4j required approximately 10 
million cycles to complete, which took approximately 48 hours on two, dual-core 2.0 GHz AMD Opteron processor. 
Clearly a well-resolved 3D disk simulation is out of the question. One way we could get around this would be to develop 
a hybrid explicit-implicit scheme, where an implicit step i s used to either solve t he radiation source terms (while using 
the current HRSC method to handle the transport) (c.f. iTurner fc Stondl2001l ) or the full radiation equations (while 
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still solving the MHD equations explicitly). Unfortunately, implementing implicit solvers in large, multidimensional 
simulations can be computationally challenging, as it involves the inversion of a large, sparsely populated matrix. It 
will take some work to determine the best way to proceed. 
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Figure 6. Luminosity I = I//^Edd as a function of mass accretion rate m = M/Mg^jj for our Bondi simulations. 
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